### Repeated Cross Section

### Implementing Method for RCS (Lebo and Weber 2015)

rm(list = ls())

pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", "fracdiff",
         "fractal",
         "lme4", "ArfimaMLM","gam",
         "forecast",
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)


# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

collapse1 <- summaryBy(aqio + aqio.lagged + aqio.lagged2 ~ day,data=p2, FUN = mean)

sd(collapse1$aqio.mean)

### checking arima first ###
# for aqi
acf(collapse1$aqio.mean) # The ACF shows that AR1 and AR2 MIGHT exist.
auto.arima(collapse1$aqio.mean) # yet, using auto.arima, it yields that p = 0, i = 1, and d = 0 (0,1,0) minimizes the AICc. 
summary(Arima(collapse1$aqio.mean, order = c(0, 1, 0))) 

auto.arima(collapse1$c_sat.mean) # 0, 0, 0
auto.arima(collapse1$l_sat.mean) # 2, 0, 2
auto.arima(collapse1$l_watch.mean) # 1, 0, 0
auto.arima(collapse1$econ.mean) # 0, 0, 0
auto.arima(collapse1$anti_west.mean) # 0, 0, 0

View(collapse1)

# effect of AQI on l_sat
a1.1 <- arfimaMLM(l_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = l_sat.mean ~ aqi.mean,
                  data = p2pp, timevar = "day")
summary(a1.1) # robust
head(a1.1$data.fd, n = 25) 

dat1 <- arfimaPrep(data = p2, timevar="day"
                   , varlist.mean = c("l_sat", "aqi","hukou","insider","soe", "educ", "kids", "married", "income",
                                  "gender", "ccp", "parade", "national", "age")
                   , varlist.fd = c("l_sat", "aqi", "parade", "national")
                   , varlist.xdif = c("hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "age"), varlist.ydif = "l_sat",  ecmformula = l_sat.mean ~ aqi.mean)
dat1 <- as.data.frame(dat1$data.merged)


fit1 <- lmer(l_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
+ educ.xdif + kids.xdif + 
  married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm, data = dat1)
summary(fit1)

fit11 <- lmer(l_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day), data = dat1)
summary(fit11)

# effect of AQI on c_sat
a1.2 <- arfimaMLM(c_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = c_sat.mean ~ aqi.mean,
                  data = p2, timevar = "day")
summary(a1.2) # robust

dat2 <- arfimaPrep(data = p2, timevar="day"
                   , varlist.mean = c("c_sat", "aqi","hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "parade", "national", "age")
                   , varlist.fd = c("c_sat", "aqi", "parade", "national")
                   , varlist.xdif = c("hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "age"), varlist.ydif = "c_sat",  ecmformula = c_sat.mean ~ aqi.mean)
dat2 <- as.data.frame(dat2$data.merged)


fit2 <- lmer(c_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm, data = dat2)
summary(fit2)

fit22 <- lmer(c_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day), data = dat2)
summary(fit22)


# effect of AQI on l_watch
a1.3 <- arfimaMLM(l_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = l_watch.mean ~ aqi.mean,
                  data = p2, timevar = "day")
summary(a1.3) # robust

dat3 <- arfimaPrep(data = p2, timevar="day"
                   , varlist.mean = c("l_watch", "aqi","hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "parade", "national", "age")
                   , varlist.fd = c("l_watch", "aqi", "parade", "national")
                   , varlist.xdif = c("hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "age"), varlist.ydif = "l_watch",  ecmformula = l_watch.mean ~ aqi.mean)
dat3 <- as.data.frame(dat3$data.merged)


fit3 <- lmer(l_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm, data = dat3)
summary(fit3)

fit33 <- lmer(l_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day), data = dat3)
summary(fit33)

# effect of AQI on c_watch
a1.4 <- arfimaMLM(c_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = c_watch.mean ~ aqi.mean,
                  data = p2, timevar = "day")
summary(a1.4) # robust

dat4 <- arfimaPrep(data = p2, timevar="day"
                   , varlist.mean = c("c_watch", "aqi","hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "parade", "national", "age")
                   , varlist.fd = c("c_watch", "aqi", "parade", "national")
                   , varlist.xdif = c("hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "age"), varlist.ydif = "c_watch",  ecmformula = c_watch.mean ~ aqi.mean)
dat4 <- as.data.frame(dat4$data.merged)


fit4 <- lmer(c_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm, data = dat4)
summary(fit4)

fit44 <- lmer(c_watch.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day), data = dat4)
summary(fit44)


# effect of AQI on anti_west
a1.6 <- arfimaMLM(anti_west.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = anti_west.mean ~ aqi.mean,
                  data = p2, timevar = "day")
summary(a1.6) # not robust

dat6 <- arfimaPrep(data = p2, timevar="day"
                   , varlist.mean = c("anti_west", "aqi","hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "parade", "national", "age")
                   , varlist.fd = c("anti_west", "aqi", "parade", "national")
                   , varlist.xdif = c("hukou","insider","soe", "educ", "kids", "married", "income",
                                      "gender", "ccp", "age"), varlist.ydif = "anti_west",  ecmformula = anti_west.mean ~ aqi.mean)
dat6 <- as.data.frame(dat6$data.merged)


fit6 <- lmer(anti_west.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm, data = dat6)
summary(fit6)

fit66 <- lmer(anti_west.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
             + educ.xdif + kids.xdif + 
               married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day), data = dat6)
summary(fit66)

# effect of AQI on life_sat
a1.7 <- arfimaMLM(life_sat.ydif ~ aqi.fd + hukou.xdif + insider.xdif + soe.xdif
                  + educ.xdif + kids.xdif + 
                    married.xdif + income.xdif + gender.xdif + ccp.xdif + parade.fd + national.fd + age.xdif + (1|day) + ecm,
                  ecmformula = life_sat.mean ~ aqi.mean,
                  data = p2, timevar = "day")
summary(a1.7) # robust

cov.labs <- c("AQI","Hukou Status","Regime Insider","SOE Employee", "Education", "Children", "Married", "Income",
                 "gender", "CCP", "Parade Period", "Holiday", "Age")
column.labels <- c("Local Government", "Central Government", "Local Oversight", "Central Oversight",  "Anti-West")

## TABLE A6
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/arfima.tex")
stargazer(fit11, fit22, fit33, fit44,fit66, font.size = "large", digits = 4, 
          column.labels = column.labels, 
          covariate.labels = cov.labs,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          star.cutoffs = c(.10, .05, .01), 
          header = F, float = F)
sink()